{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These examples are based on **Chapter 15** of _Introduction to Econometrics_ by Jeffrey Wooldridge and demonstrate the basic use of the IV estimators (primarily IV2SLS -- the two-stage least squares estimator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wages of Married Women\n",
    "\n",
    "This first example examines the effect of education on the wages of women.  Education is a classic endogenous variable since it has signaling value beyond the actual direct effect (among other reasons).\n",
    "\n",
    "This first block imports the data and uses the `DESCR` attribute to describe the dataset.  `add_constant` from statsmodels is used to add a variable named `const` to the DataFrame.  Observations with missing values are dropped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from statsmodels.api import add_constant\n",
    "from linearmodels.datasets import mroz\n",
    "\n",
    "print(mroz.DESCR)\n",
    "data = mroz.load()\n",
    "data = data.dropna()\n",
    "data = add_constant(data, has_constant=\"add\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since OLS is a special case of 2SLS, IV2SLS can be used to estimate a model using OLS by setting `endog` and `instruments` variables to `None`.\n",
    "\n",
    "This first regression uses OLS to estimate the effect of education on the log of wage. It indicates that 1 year of education increases wage by 10%.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from linearmodels.iv import IV2SLS\n",
    "\n",
    "res_ols = IV2SLS(np.log(data.wage), data[[\"const\", \"educ\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two-stage least squares is estimated \"as-if\" two regressions are run.  Here the first stage regression where OLS is used to fit the value on the instrument (in this case the education of the subject\"s father). The fitted value is saved for use later.\n",
    "\n",
    "This first stage regression indicates that there is a strong relationship and the first stage easily passes the rule of thumb test where the F-statistic is at least 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_first = IV2SLS(data.educ, data[[\"const\", \"fatheduc\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res_first)\n",
    "data[\"educ_hat\"] = data.educ - res_first.resids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second stage uses the instrument to fit the model.  This model indicates a much lower effect of education. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_second = IV2SLS(np.log(data.wage), data[[\"const\"]], data.educ, data.fatheduc).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res_second)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitted value can be used with OLS to estimate the same parameters.  Note that the other values reported such as t-statistics or $R^2$ are not correct.\n",
    "\n",
    "The `compare` function is used to compare the different models. The 2SLS coefficient on education and the direct coefficient on educ_hat are identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.iv import compare\n",
    "\n",
    "res_direct = IV2SLS(np.log(data.wage), data[[\"const\", \"educ_hat\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(compare({\"OLS\": res_ols, \"2SLS\": res_second, \"Direct\": res_direct}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wages of Men\n",
    "This next example examines the returns to education for men and uses the number of siblings as an instrument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.datasets import wage\n",
    "\n",
    "men = wage.load()\n",
    "print(wage.DESCR)\n",
    "men = men[[\"educ\", \"wage\", \"sibs\", \"exper\"]]\n",
    "men = add_constant(men)\n",
    "men = men.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OLS estimates indicate a small effect of education for men in this sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_ols = IV2SLS(np.log(men.wage), men[[\"const\", \"educ\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first stage regressed the endogenous variable on the instrument.  There is a strong, negative relationship here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_first = IV2SLS(men.educ, men[[\"const\", \"sibs\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res_first.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second stage indicates a much strong relationship than the OLS estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(np.log(men.wage), men.const, men.educ, men.sibs).fit(cov_type=\"unadjusted\")\n",
    "print(res.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weighted IV\n",
    "\n",
    "All IV estimator support weighted which extends the concept of WLS in an OLS model to IV estimation.  The weights are applied to the dependent variables, the matrix of regressors (endogenous and exogenous) and the matrix of instrument (exogenous and instruments) which allows for GLS-type estimation.  In particular, if the variance of model residuals was $\\sigma^2_i$, then setting $w_i = 1 / \\sigma^2_i$ would produce GLS estimates. \n",
    "\n",
    "This example shows how a feasible GLS-type of estimator could be used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res2 = res.resids ** 2\n",
    "fgls_mod = IV2SLS(np.log(res2), men[[\"const\", \"sibs\", \"exper\"]], None, None)\n",
    "fgls_res = fgls_mod.fit()\n",
    "sigma2_hat = np.exp(np.log(res2) - fgls_res.resids)\n",
    "print(fgls_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance of the squared residuals is not particularly well explained by the data, and so the GLS-type estimates and the usual IV estimates don\"t differ by much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_gls = IV2SLS(\n",
    "    np.log(men.wage), men.const, men.educ, men.sibs, weights=1 / sigma2_hat\n",
    ").fit(cov_type=\"unadjusted\")\n",
    "print(res_gls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smoking and birth weight\n",
    "This example examines the effect of smoking on the birth weight of babies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.datasets import birthweight\n",
    "\n",
    "data = birthweight.load()\n",
    "print(birthweight.DESCR)\n",
    "data = add_constant(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first stage regresses the number of packs smoked on the cigarette price.  There is a very weak relationship -- so weak that this is an example of a _weak instrument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(data.packs, data[[\"const\", \"cigprice\"]], None, None).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the weak relationship between the price and the number of pack smoked, the 2SLS can be estimated, although substantial caution is warranted to interpret the results. Note the very negative $R^2$ (-150%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(np.log(data.bwght), data.const, data.packs, data.cigprice).fit(\n",
    "    cov_type=\"unadjusted\"\n",
    ")\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proximity and education\n",
    "\n",
    "This next example uses a well-known dataset that uses proximity to a 4 year college as an instrument for education."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.datasets import card\n",
    "\n",
    "data = card.load()\n",
    "print(card.DESCR)\n",
    "data = add_constant(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lists are used to hold the groups of variable in this large model and missing values are dropped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dep = [\"wage\"]\n",
    "endog = [\"educ\"]\n",
    "exog = [\n",
    "    \"const\",\n",
    "    \"exper\",\n",
    "    \"expersq\",\n",
    "    \"black\",\n",
    "    \"smsa\",\n",
    "    \"south\",\n",
    "    \"smsa66\",\n",
    "    \"reg662\",\n",
    "    \"reg663\",\n",
    "    \"reg664\",\n",
    "    \"reg665\",\n",
    "    \"reg666\",\n",
    "    \"reg667\",\n",
    "    \"reg668\",\n",
    "    \"reg669\",\n",
    "]\n",
    "instr = [\"nearc4\"]\n",
    "data = data[dep + exog + endog + instr].dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first stage estimate shows a very large F-statistic.  Note that when there are many exogenous variables the results cannot be directly interpreted.  It is better to use the `first_stage` information from a 2SLS estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(data.educ, data[instr + exog], None, None).fit()\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The OLS estimate indicates an increase of 7% for a year of education."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(np.log(data.wage), data[exog + endog], None, None).fit()\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 2SLS estimate of the effect is nearly double.  However, there is some reason to be concerned about the strength of the instrument despite the F-statistic in the first stage regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_2sls = IV2SLS(np.log(data.wage), data[exog], data[endog], data[instr]).fit()\n",
    "print(res_2sls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The property `first_stage` can be used to show a large set of first stage diagnostics. These results show a much lower _partial_ $R^2$ that has measures the unique effect of the instrument on the endogenous controlling for the exogenous regressors.  This is much smaller than the naive first stage $R^2$ of 47%.  The partial F-statistic is also much smaller, although it technically over the rule-of-thumb of 10 for a single instrument.\n",
    "\n",
    "The instrument is a dummy variable and being close to a 4 year college is only worth 0.3 years of education on average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_2sls.first_stage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Formula interface\n",
    "\n",
    "This model was large and so it might be simpler to use a formula. While formulas are discussed in [detail in another notebook](using-formulas.ipynb), they use the formula language of patsy with an augmentation to specify the endogenous and instrumental variables. The generic form is\n",
    "\n",
    "```\n",
    "dependent ~ exog + [endog ~ instr]\n",
    "```\n",
    "\n",
    "where each block can contain multiple variables.\n",
    "\n",
    "Here the model is compared to the direct parameterization using DataFrames by differencing the $R^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "formula = (\n",
    "    \"np.log(wage) ~ 1 + exper + expersq + black + smsa + south + smsa66 + reg662 + reg663 + reg664 + \"\n",
    "    \"reg665 + reg666 + reg667 + reg668 + reg669 + [educ ~ nearc4]\"\n",
    ")\n",
    "mod = IV2SLS.from_formula(formula, data)\n",
    "res_formula = mod.fit(cov_type=\"unadjusted\")\n",
    "print(res_formula.rsquared - res_2sls.rsquared)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Categorical Variables\n",
    "pandas ``categorical``s are automatically treated as factors and expanded to dummies.  The first is always dropped. This next block constructs a categorical from the region dummies and then uses it instead of the individual dummies.  The model is identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"reg\"] = \"661\"  # The default region, which was omitted\n",
    "for i in range(2, 10):\n",
    "    region = \"reg66\" + str(i)\n",
    "    data.loc[data[region] == 1, \"reg\"] = region[3:]\n",
    "data[\"reg\"] = data[\"reg\"].astype(\"category\")\n",
    "data.describe()\n",
    "res_cat = IV2SLS(\n",
    "    np.log(data.wage),\n",
    "    data[[\"const\", \"exper\", \"expersq\", \"black\", \"smsa\", \"south\", \"smsa66\", \"reg\"]],\n",
    "    data.educ,\n",
    "    data.nearc4,\n",
    ").fit()\n",
    "print(res_cat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Post-estimation diagnostics\n",
    "\n",
    "Common post-estimation diagnostics are to test the assumption of endogeneity and to examine if instruments are valid (when there are more instruments than endogenous variables)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = mroz.load()\n",
    "data = data.dropna()\n",
    "data = add_constant(data, has_constant=\"add\")\n",
    "data[\"lnwage\"] = np.log(data.wage)\n",
    "dep = \"lnwage\"\n",
    "exog = [\"const\", \"exper\", \"expersq\"]\n",
    "endog = [\"educ\"]\n",
    "instr = [\"fatheduc\", \"motheduc\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to fit the model using 2SLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = IV2SLS(data[dep], data[exog], data[endog], data[instr]).fit(cov_type=\"unadjusted\")\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wooldridge\"s regression test of exogeneity uses regression residual where the endogenous variables are regressed on the exogenous and the instrument to test for endogeneity. IF the endogenous variable is actually exogenous these residuals should not be correlated with the variable of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.wooldridge_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This test can be easily implemented using two regression.  The first one constructs the residuals and the second re-fits the model using 2SLS but including the residuals. \n",
    "\n",
    "Note that the p-value of the t-state on `residuals` is the same as the P-value of the previous test -- this is not an accident."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = IV2SLS(data[endog], data[exog + instr], None, None).fit().resids\n",
    "import pandas as pd\n",
    "\n",
    "res_direct = IV2SLS(\n",
    "    data[dep], pd.concat([v, data[exog]], 1), data[endog], data[instr]\n",
    ").fit(cov_type=\"unadjusted\")\n",
    "print(res_direct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this regression has two instrument it is possible to test for overidentification.  Wooldridge\"s overidentification test uses a regression to test whether the 2SLS residuals are uncorrelated with the instruments, which should be the case if the model is correct and the instruments are not needed in the original model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.wooldridge_overid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A naive version of this test can be directly implemented. This direct implementation is different from the formal test but would be consistent if the model was overidentified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = res.resids\n",
    "res = IV2SLS(u, data[[\"exper\", \"expersq\"] + instr], None, None).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test is $n\\times R^2$, and has the same $\\chi^2_1$ distribution. The test statistic is slightly smaller but the conclusions are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.nobs * res.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Husband\"s education can be used as an additional instrument, and its validity tested."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instr = [\"fatheduc\", \"motheduc\", \"huseduc\"]\n",
    "res = IV2SLS(data[dep], data[exog], data[endog], data[instr]).fit(cov_type=\"unadjusted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing overidentification does not indicate any difference from the previous result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.wooldridge_overid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directly testing using two regression would reach the same conclusion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = res.resids\n",
    "res = IV2SLS(u, data[[\"exper\", \"expersq\"] + instr], None, None).fit()\n",
    "res.nobs * res.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Panel IV \n",
    "Instrumental variable regression can also be used with panel data. This example makes use of first differences to eliminate a year-specific effect and then uses "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.datasets import jobtraining\n",
    "\n",
    "data = jobtraining.load()\n",
    "print(jobtraining.DESCR)\n",
    "data.head()\n",
    "data = data.where(data.year.isin((1987, 1988)))\n",
    "data = data.dropna(how=\"all\", axis=0).sort_values([\"fcode\", \"year\"])\n",
    "print(data.describe())\n",
    "data = data.set_index(\"fcode\")\n",
    "data = data[[\"year\", \"hrsemp\", \"grant\", \"scrap\", \"lscrap\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltas = data.loc[data.year == 1988] - data.loc[data.year == 1987]\n",
    "deltas = add_constant(deltas, has_constant=\"add\")\n",
    "deltas = deltas.dropna()\n",
    "print(deltas.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first stage indicates a relatively strong relationship between grant and the number of hours employed.  Note that grant is a dummy and so the coefficient is just the difference in means. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = IV2SLS(deltas.hrsemp, deltas[[\"const\", \"grant\"]], None, None)\n",
    "print(mod.fit(cov_type=\"unadjusted\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here a formula is used to specify the model since it is cleaner.  Note that the `[]` contains the endogenous variables and instruments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = IV2SLS.from_formula(\"lscrap ~ 1 + [hrsemp ~ grant]\", deltas)\n",
    "res_iv = mod.fit(cov_type=\"unadjusted\")\n",
    "print(res_iv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 2SLS estimate is nearly twice as large as the OLS estimate and is slightly more significant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_ols = IV2SLS.from_formula(\"lscrap ~ 1 + hrsemp\", deltas).fit(cov_type=\"unadjusted\")\n",
    "print(compare({\"Panel OLS\": res_ols, \"Panel IV\": res_iv}))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
